!> \brief Contains tools subroutines for the main algorithm.
module tools

   implicit none

contains

   !============================================================================

   !> \brief Generates a random initial population.
   subroutine get_initial_population(g, nu, np, xmin, xmax, pop)
      implicit none
      integer, intent(out) :: g            !< generation
      integer, intent(in)  :: nu           !< number of unknowns
      integer, intent(in)  :: np           !< population size
      real(8), intent(in)  :: xmin(nu)     !< low boundary constraints
      real(8), intent(in)  :: xmax(nu)     !< high boundary constraints
      real(8), intent(out) :: pop(np,nu)   !< population

      integer :: j
      integer :: ind
      real(8) :: rnd

      g = 1

      do ind = 1, np

         do j = 1, nu
            call random_number(rnd)
            pop(ind,j) = xmin(j) + rnd*(xmax(j) - xmin(j))
         end do

      end do

   end subroutine get_initial_population


   !> \brief Calculates a random individual
   subroutine get_random_individual(nu, xmin, xmax, x)
      implicit none
      integer, intent(in) :: nu           !< number of unknowns
      real(8), intent(in) :: xmin(nu)     !< low boundary constraints
      real(8), intent(in) :: xmax(nu)     !< high boundary constraints
      real(8), intent(out) :: x(nu)       !< random individual

      integer :: j   ! dummy index
      real(8) :: rnd ! random number

      do j = 1, nu

         call random_number(rnd)

         x(j) = xmin(j) + rnd*(xmax(j) - xmin(j))

      end do

   end subroutine get_random_individual


   !============================================================================

   !> \brief Selects three distinct individuals of the population.
   subroutine select_individuals(np, ind, r)
      implicit none
      integer, intent(in) :: np     !< population size
      integer, intent(in) :: ind    !< number of the individual
      integer, intent(out) :: r(3)  !< indices of selected individuals

      real(8) :: rnd

      r = ind

      do while (r(1) == ind)
         call random_number(rnd)
         r(1) = int(rnd*np) + 1
      end do

      do while (r(2) == r(1) .or. r(2) == ind)
         call random_number(rnd)
         r(2) = int(rnd*np) + 1
      end do

      do while (r(3) == r(1) .or. r(3) == r(2) .or. r(3) == ind)
         call random_number(rnd)
         r(3) = int(rnd*np) + 1
      end do

   end subroutine select_individuals

   !============================================================================

   !> \brief Generates a trial individual from randomly selected individuals
   !! of the population.
   subroutine get_trial_individual(ind, nu, np, dif, crs, r, pop, x)
      implicit none
      integer, intent(in) :: ind          !< number of the individual
      integer, intent(in) :: nu           !< number of unknowns
      integer, intent(in) :: np           !< population size
      integer, intent(in) :: r(3)         !< indices of selected individuals
      real(8), intent(in) :: dif          !< differentiation constant
      real(8), intent(in) :: crs          !< crossover constant
      real(8), intent(out) :: x(nu)       !< trial individual
      real(8), intent(in) :: pop(np,nu)   !< population

      integer :: j
      integer :: irnd
      real(8) :: rnd

      call random_number(rnd)

      irnd = int(rnd*nu) + 1

      do j = 1, nu

         call random_number(rnd)

         if (rnd < crs .or. irnd == j) then
            x(j) = pop(r(3),j) + dif*(pop(r(1),j) - pop(r(2),j))
         else
            x(j) = pop(ind,j)
         end if

      end do

   end subroutine get_trial_individual

   !============================================================================

   !> \brief Verifying the boundary constraints
   subroutine get_constraints_check(nu, xmin, xmax, x)
      implicit none
      integer, intent(in) :: nu           !< number of unknowns
      real(8), intent(inout) :: x(nu)     !< trial individual
      real(8), intent(in) :: xmin(nu)     !< low boundary constraints
      real(8), intent(in) :: xmax(nu)     !< high boundary constraints

      integer :: j
      real(8) :: rnd

      do j = 1, nu

         if (x(j) < xmin(j) .or. x(j) > xmax(j)) then

            call random_number(rnd)
            x(j) = xmin(j) + rnd*(xmax(j) - xmin(j))

         end if

      end do

   end subroutine get_constraints_check

   !============================================================================

   !> \brief Get the fitness of each individual
   subroutine get_fitness(folderout, sname, ffit, ind, nu, x, xname, xfit, estatus) ! last parameter is output
      implicit none
      character(200), intent(in) :: folderout   !< folder the for output files
      character(200), intent(in) :: sname       !< simulations name
      character(200), intent(in) :: ffit        !< name of executable for fitness calculation
      integer, intent(in) :: ind    !< individual number
      integer, intent(in) :: nu     !< dimension
      real(8), intent(in) :: x(nu)  !< parameters for optimization
      character(10), intent(in) :: xname(nu) !< names of the parameters
      real(8), intent(out) :: xfit  !< fitness
      integer, intent(out) :: estatus !< exit status ( 0 = success, 1 = failure, 2 = generate another individual)

      character(3) :: char3
      character(200) :: arqpar   !< name of file for the parameters
      character(200) :: arqfit   !< name of file for the fitness
      integer :: i

      call convert_int_to_char3(ind, char3)

      arqpar = trim(folderout) // "parameters" // char3 // ".txt"
      arqfit = trim(folderout) // "fitness"    // char3 // ".txt"

      open(27, file = trim(arqpar))

      if (len(trim(arqfit)) <= 23) then
         write(27,"(a23, a)") "'" // trim(arqfit) // "'", " = arqfit: name of the file for fitness"
      else
         write(27,"(a  , a)") "'" // trim(arqfit) // "'", " = arqfit: name of the file for fitness"
      end if

      if (len(trim(sname)) <= 23) then
         write(27,"(a23, a)") trim(sname), " = sname: simulation name"
      else
         write(27,"(a  , a)") trim(sname), " = sname: simulation name"
      end if

      write(27,"(i23, a)") ind, " = ind:   individual number"
      write(27,"(i23, a)")  nu, " = nu:    number of unknowns"

      write(27,*)
      write(27,*) "=========== VARIABLES OF OPTIMIZATION ====================="
      write(27,*) "                     X                       Variable name"

      do i = 1, nu
         write(27,"(1pe23.15, a36)") x(i), trim(xname(i))
      end do

      write(27,*)
      close(27)

      call system("./" // trim(ffit) // " < " // trim(arqpar) // " > /dev/null")


      open(26, file = trim(arqfit))
      read(26,*) xfit
      read(26,*) estatus
      close(26)

   end subroutine get_fitness

   !============================================================================

   !> \brief Converts an integer to a string of length 3.
   subroutine convert_int_to_char3(nmbr, char3) ! last parameter is output
      implicit none
      character(3), intent(out) :: char3
      integer, intent(in) :: nmbr

      if (nmbr < 10) then
         write(char3,"('00',i1)") nmbr
      else if (nmbr < 100) then
         write(char3,"('0',i2)") nmbr
      else
         write(char3,"(i3)") nmbr
      end if

   end subroutine convert_int_to_char3

   !============================================================================


   !> \brief Calculates the convergence measure based on several methods
   !! kcm = 1: P-measure with dimensional variables
   !! kcm = 2: P-measure with dimensionless variables
   subroutine get_convergence(kcm, nu, np, xmin, xmax, pop, fit, cm)
      implicit none
      integer, intent(in)  :: kcm        !< Kind of convergence measure
      integer, intent(in)  :: nu         !< Dimension of the problem
      integer, intent(in)  :: np         !< Size of the population
      real(8), intent(in)  :: xmin(nu)   !< low boundary constraints
      real(8), intent(in)  :: xmax(nu)   !< high boundary constraints
      real(8), intent(in)  :: pop(np,nu) !< Population
      real(8), intent(in)  :: fit(np)    !< Fitness function
      real(8), intent(out) :: cm         !< Convergence measure


      select case (kcm)

         case (1) ! P-measure with dimensional variables

            call get_p_measure( 1, nu, np, xmin, xmax, pop, cm)

         case (2) ! P-measure with dimensionless variables

            call get_p_measure( 2, nu, np, xmin, xmax, pop, cm)

         case default

            write(*,*) "ERROR: unknown kind of convergence measure."

            stop

      end select

   end subroutine get_convergence

   !============================================================================

   !> \brief Calculates the P-measure, i. e. the population
   !! convergence measure of FEOKTISTOV "Differential evolution" 2006
   subroutine get_p_measure(kpm, nu, np, xmin, xmax, pop, pm)
      implicit none
      integer, intent(in)  :: kpm        !< Kind of P-measure (1=dimensional, 2=dimensionless)
      integer, intent(in)  :: nu         !< Dimension of the problem
      integer, intent(in)  :: np         !< Size of the population
      real(8), intent(in)  :: xmin(nu)   !< low boundary constraints
      real(8), intent(in)  :: xmax(nu)   !< high boundary constraints
      real(8), intent(in)  :: pop(np,nu) !< Population
      real(8), intent(out) :: pm         !< P-measure

      ! Inner variables

      integer :: j ! Dummy index

      real(8) :: raux ! Auxiliary variable

      real(8), dimension(nu)    ::   pb ! Mean individual
      real(8), dimension(np,nu) :: apop ! Auxiliary population


      if ( kpm == 1 ) then

         ! Dimensional population

         apop = pop

      else

         ! Calculation of the dimensionless population

         do j = 1, np

            apop(j,:) = ( pop(j,:) - xmin ) / ( xmax - xmin )

         end do

      end if

      ! Calculating the mean individual

      pb = 0.d0

      do j = 1, np

         pb = pb + apop(j,:)

      end do

      pb = pb / dble(np)


      ! Calculating the P-measure

      pm = 0.d0

      do j = 1, np

         raux = norm2( apop(j,:) - pb )

         if ( raux > pm ) pm = raux

      end do

   end subroutine get_p_measure

   !============================================================================

end module tools
